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In several environmental applications data are functions of time, essentially continuous, 
observed and recorded discretely, and spatially correlated. Most of the methods for an- 
alyzing such data are extensions of spatial statistical tools which deal with spatially 
dependent functional data. In such framework, this paper introduces a new clustering 
method. The main features are that it finds groups of functions that are similar to each 
other in terms of their spatial functional variability and that it locates a set of centers 
which summarize the spatial functional variability of each cluster. The method opti- 
mizes, through an iterative algorithm, a best fit criterion between the partition of the 
curves and the representative element of the clusters, assumed to be a variogram func- 
tion. The performance of the proposed clustering method was evaluated by studying 
the results obtained through the application on simulated and real datasets. 
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Introduction 



Spatial interdependence of phenomena is a common feature of many environmental applications such 
as oceanography, geochemistry, geometallurgy, geography, forestry, environmental control, landscape 
ecology, soil science, and agriculture. For instance, in daily patterns of geophysical and environmental 
phenomena where data (from temperature to sound) are instantaneously recorded over large areas, ex- 
planatory variables are functions of time, essentially continuous, observed and recorded discretely, and 
spatially correlated. 

In the last years, the analysis of such data has been performed by Spatial Functional Data Analysis 
(SFDA) (Delicado et al. (2010)), a new branch of Functional Data Analysis (Ramsay, Silverman (2005)). 
Most of the contributions in this framework are extensions of spatial statistical tools for functional data. 

This paper focuses on clustering spatially related curves. 

To the authors knowledge, existent clustering strategies for spatially dependent functional data are 
very limited. The approaches refer to the following main methods: hierarchical, dynamic, clusterwise 
and model-based. The hierarchical group of methods, (Giraldo et al. (2009)) is based on spatial weighted 
dissimilarity measures between curves. These are extensions to the functional framework of the ap- 
proaches proposed for geostatistical data, where the norm between curves is replaced by a weighted 
norm among the geo-referenced functions. In particular, two alternatives are proposed for univariate and 
multivariate context, respectively. In the univariate framework, the weights correspond to the variogram 
values computed for the distance between the sites. In the multivariate framework, a dimensionality re- 
duction is performed using a Principal Component Analysis technique for functional data (Dauxois et 
al. (1982)) with the variogram values, computed on the first principal component, used as weights. The 
main characteristic of these approaches is in considering the spatial dependence among different kinds of 
functional data and in defining spatially weighted distances measures. 

Alternatively to these approaches, with the aim of obtaining a partition of spatial functional data and 
a suitable representation for each cluster, the same authors proposed dynamic (Romano et al. (2010)) 
and clusterwise methods (Romano, Verde (2009)). The first, aims at classifying spatially dependent 
functional data and achieving a kriging spatio-functional model prototype for each cluster by minimizing 
the spatial variability measure among the curves in each cluster. 
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In the ordinary kriging for functional data, the problem is to obtain an estimated curve in an unsam- 
pled location. This proposed method gets not only a prediction of the curve but also a best representative 
location. In this sense, the location is a parameter to estimate and the objective function may have sev- 
eral local minima corresponding to different local kriging. The method proposes to solve this problem 
by evaluating local kriging on unsampled locations of a regular spatial grid in order to obtain the best 
representative predictor for each cluster. This approach is based on the definition of a grid of sites in 
order to obtain the best representative function. In a different manner and for several functional data, the 
clusterwise linear regression approach attempts to discover spatial functional linear regression models 
with two functional predictors, an interaction term, and spatially correlated residuals. This approach can 
establish a spatial organization in relation to the interaction among different functional data. The algo- 
rithm is a k-means clustering with a criterion based on the minimization of the squared residuals instead 
of the classical within cluster dispersion. 

A further approach is a model-based method for clustering multiple curves or functionals under spa- 
tial dependence specified by a set of unknown parameters (Jiang, Serban (2010)). The functionals are 
decomposed using a semi-parametric model, with fixed and random effects. The fixed effects account 
for the large-scale clustering association and the random effects account for the small scale spatial de- 
pendence variability. Although the clustering algorithm is one of the first endeavors in handling densely 
sampled space domains using rigorous statistical modeling, it presents several computational difficulties 
in applying the estimation algorithm to a large number of spatial units. 

The method proposed in this paper, belongs to the dynamic clustering approaches (Diday (1971)). 

The current interest is motivated by a wide number of environmental applications where understanding 

the spatial relation among curves in an area is an important source of information for making a prediction 

regarding an unknown point of the space. The main idea is to provide a summary of the set of curves 

spatially correlated by a prototype-based clustering approach. With this aim the proposed method uses a 

Dynamic Clustering approach to optimize a best fit criterion between the partition and the representative 

element of the clusters, assumed to be a variogram function. [] According to this procedure, clusters 

are groups of functions that are similar to each other in terms of their spatial functional variability. The 

central issue in the procedure consists in taking into account the spatial dependence of georeferenced 
'A preliminary version of this paper appears in (Romano et al. (2010)) 
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functional data. For most environmental applications, the spatial process is considered to be stationary 
and isotropic, and a wide area of the space is modeled with a single variogram model. In practice, 
however, many spatial functional data cannot be modeled accurately with the same variogram model. 
Recognizing this, the scope is to propose a clustering method that clusters the geo-referenced curves into 
groups and associates a variogram function to each of them. 

The rest of this paper is organized as follows. Section 1 introduces the concept of spatial functional 
data and the measures for studying their spatial relation. Section 2 shows the proposed method. Section 
3 illustrates the method on synthetic and real datasets. 

1 Spatial variability measure for geostatistical functional data 

Spatially dependent functional data may be defined as the data for which the measurements on each 
observation that is a curve are part of a single underlying continuous spatial functional process defined as 

{Xs ■■ seDC R d } (1) 

where s is a generic data location in the d— dimensional Euclidean space (d is usually equal to 2), the set 
D C R d can be fixed or random, and Xs are functional random variables, defined as random elements 
taking values in an infinite dimensional space. The nature of the set D allows the classification of Spatial 
Functional Data. Following (Delicado et al. (2010)) these can be distinguished in geostatistical functional 
data, functional marked point patterns and functional areal data. 

The paper focuses on geostatistical functional data, where samples of functions are observed in dif- 
ferent sites of a region (spatially correlated functional data). 

Let {xs(t) ■ t G T, s G D C R d ) be a random field where the set D C R d is a fixed subset of R d 
with positive volume. Xs is a functional variable defined on some compact set T of R for any s G D. 

It is assumed to observe a sample of curves (Xsi(t), ■ ■ ■ , Xs x (t), ■ ■ ■ , Xs„(t)) for t G T where Sj is a 
generic data location in the rf-dimensional Euclidean space. 

For each t, the random process is assumed to be second order stationary and isotropic: that is, the 
mean and variance functions are constant and the covariance depends only on the distance between sam- 
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pling sites. Formally: E(x 8 (t)) = m(t), for all t G T, s G D, V(x«(f)) = <r 2 (*), for all t G T, s G A 
and Cou(x Si (i),x s .(t)) = C(h,t) where = ||sj — sj|| and all s i: Sj G -D 

This implies that a variogram function for functional data "f(h, t) exists, also called trace-variogram 
function (Giraldo et al. (2009)), such that 



7 (M) = J Si sM = l^iXsM ~ XsM = \^ [XsM ~ XsMT ( 2 ) 

where h — ||sj — Sj || and all Sj, s^- G £>. 

By using Fubini's theorem, the previous becomes j(h) = § T ^ SiS .(t)dt for \\si — Sj\\ = h. This 
variogram function can be estimated by the classical method of the moments by means of: 



1 V 71 i,jeN(h) JT 



(3) 



where iV(/i) = {(s,; Sj) : ||s, — Sj\\ = h} for regular spaced data and \N(h)\ is the number of distinct 
elements in N(h). 

When data are irregularly spaced, N(h) = {(s^ Sj) : \\si — Sj\\ G (h — e,h + e)} with e > being a 
small value. 

The estimation of the empirical variogram for functional data using ([3]) involves the computation of 
integrals that can be simplified by considering that the functions are expanded in terms of some basis 
functions 

z 

Xs i (t) = ^2a a B l (t) = a i T B(t), i = l,...,n (4) 



i=i 



where a* is the vector of the basis coefficients for the % Si , then the coefficients of the curves can be 
consequently organized in a matrix as follows: 



^ ai i ai 2 • • • ai z ^ 
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\ 0>n,l a n,2 ■ ■ ■ &2,Z J 
\ ' ' / nxZ 

Thus, the empirical variogram function for functional data can be obtained by considering: 



/ {XsAt) - XsM 2 dt = [ (a, T B(t)-a/B(t)) 2 ^ = 
Jt Jt 

= J (a i -a J ) T B(t) 2 dt = 




= (a l -a,) T W(a i -a,) T 

where W = f T B (t) B(t) T dt is the Gram matrix that is the identity matrix for any orthonormal basis. 
For other basis as B-Spline basis function, W is computed by numerical integration. Thus the variogram 
is expressed by: 

7 ^ = 2\N(h)\ ^ ~ W ^ ~ V '' 3 1 11 Si ~ S ^ = h 

1 1 )l i,jeN(h) 

The empirical variograms cannot be computed at every lag distance h, and due to variation in the 
estimation, it is not ensured that it is a valid variogram. 

In applied geostatistics, the empirical variograms are thus approximated (by ordinary least squares 
(OLS) or weighted least squares (WLS)) by model functions, ensuring validity (Chiles, Delfiner (1999)). 
Some widely used models include: Spherical, Gaussian, exponential, or Mathern (Cressie (1993)). The 
variogram, as defined before, is used to describe the spatial variability among functional data across an 
entire spatial domain. In this case, all possible location pairs are considered. 

However, this spatial variability may be strongly influenced by an unusual or changing behavior 
within this wide area. For instance, in climatology, a sensor network is used to evaluate the temperature 
variability over an area. Some sensors could describe the characteristics of their surrounding sites with 
very different proportions, causing potentials for errors in the computation of spatial variability. 

Thus, in order to describe these spatial variability substructures, this paper introduces the concept of 
the spatial variability components with regards to a specific location by defining a centered variogram for 
functional data. 

Coherently with the above definition, given a curve Xsi(0> me centered variogram for functional data 
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can be expressed by 

7*(M) = ^E(X«(*) -*.,(*)) (5) 

for each Sj ^ Si E D. Similar to the variogram function, the centered variogram of the curve x Si (t), as a 
function of the lag h, can be estimated through the method of moments: 

'•'• W = 2W^)\ E J T (x, i it)-X, l (Ofdt (6) 

where N s *(h) C iV(/i) = Sj) : \\si - Sj\\ = h} and it is such that \N(h)\ = £\ \N Si (h)\. 
Through straightforward algebraic operations, it is possible to show that the variogram function is a 
weighted average of centered variograms: 

thus: 

1 ™ 

It is worth noting that the estimation of the centered variogram can be expressed in the same manner 
in the functional setting. 



2 Variogram-based Dynamic Clustering approach for spatially de- 
pendent functional data 

A Dynamic Clustering Algorithm (DC A) (Celeux et al. (1988)) (Diday (1971)) is an unsupervised learn- 
ing algorithm, which finds partitions a set of objects into internally dense and sparsely connected clusters. 
The main characteristic of the DCA is that it finds, simultaneously, the partition of data into a fixed num- 
ber of clusters and a set of representative syntheses, named prototypes, obtained through the optimization 
of a fitting criterion. Formally, let E be a set of n objects. The Dynamic Clustering Algorithm finds a 
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partition P* = (d, 
L* = (G 1 , . . . , G k , 
criterion: 



. , Ck, ■ ■ ■ , Ck) of E in K non empty clusters and a set of representative prototypes 
. , G K ) for each C k cluster of P so that both P* and L* optimize the following 



A(P*, L*) = Mm {A(P, L) / PeP k ,Le A K } (9) 

with P^ the set of all the ^-cluster partitions of E and the representation space of the prototypes. 
A(P, L) is a function, which measures how well the prototype Gk represents the characteristics of objects 
of the cluster and it can usually be interpreted as an heterogeneity or a dissimilarity measure of goodness 
of fit between Gk and C k . 

The definition of the algorithm is performed according to two main tasks: 

- representation function allowing to associate to each partition P E Pr of the data in K classes C k 
(k — 1, . . . , K), a set of prototype L = (G 1: . . . , G k , . . . , G K ) of the representation space 

- allocation function allowing to assign to each G k E L, a set of elements C k . 

The first choice concerns the representation structure L for the classes C\, . . . , Ck G P. 

Let {t),..., Xs n (t)} (with t E T and s G P) be the sample of spatially located functional data. 
The proposed method aims at partitioning them into clusters in order to minimize, in each cluster, the 
spatial variability. 

Following this aim, the method optimizes a best fit criterion between the centered variogram function 
7^(/i) and a theoretical variogram function ^(h) for each cluster as follows: 

HP,L) = J2 E (tfW-TJW) 2 (10) 

where 7^ is the centered variogram, which describes the spatial dependence between a curve Xs x (t) at 
the site s, and all the other curves (t) at different spatial lags ft,. This allows to evaluate the membership 
of a curve x Si (t) to the spatial variability structure of an area. 

As already mentioned, starting from a random initialization, the algorithm alternates representation 
and allocation steps until it reaches the convergence to a stationary value of the criterion A(P, L). 
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In the representation step, the theoretical variogram 'jl(h) of the set of curves Xs»(0 e Cfe> f° r eacri 
cluster C fc is estimated. This involves the computation of the empirical variogram and its model fitting 
by the Ordinary Least Square method. 

In the allocation step, the function 7^ is computed for each curve Xsi(t)- Then a curve XsAt) is 
allocated to a cluster C k by evaluating its matching with the spatial variability structure of the clusters 
according to the following rule: 

E (T? CO - 72W)% < E WW - T*'C0)V VA; ^ fc' (11) 

h<h*£{rn k ;M k ] h<h* d[m k ;Mk\ 

where: 

• Pk — |jv fc | and p fe ' = MM- are the weights computed respectively, considering, for a fixed s*, the 
number of location pairs A/"*, 1 that are separated by a distance /i in a cluster /c, and A;'. 

• m k — min fc /i^, , M k = max fc /i^ where h* k is the spatial distance at which the variogram 7^ for each 
cluster k reaches its sill. 

The problem is that for each cluster, there are several values of h* k (k — 1, . . . , K), due to the different 
spatial functional variability structures of the partition. According to the above allocation criterion, only 
one level h* is chosen such that for h > h* , there is no spatial correlation. This rule facilitates the spatial 
aggregation process leading to a tendency to form regions of spatially correlated curves. Especially, h* 
is set in the range [m k , M k \. 

The consistency between the representation of the clusters and the allocation criterion guarantees the 
convergence of the criterion to a stationary minimum value (Celeux et al. (1988)). 

In the context of the proposed method, this is verified when: 

11(h) = argmin E (tfCO - 72C0) 2 ( 12 ) 
x Si (t)ec k 

Thus, since the allocation of each curve Xs l (t) to a cluster C k is based on computing the squared 
Euclidean distance between ^(h) and ^(h), since the variogram ^(h) is the average of the functions 
7^(/i), then 7j*(/i) minimizes the spatial variability of each cluster. 
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Algorithm 1 Dynamic Clustering Algorithm for geostatistical functional data 
Initialization: 

Start from a random partition P = (Ci, . . . , . . . , CV) 
Representation step: 
for all clusters C fc do 

Compute the prototype 7j£(/i) which optimizes the best fitting criterion: 

mm £ (7?(/0-72C0) 2 

x Si (t)eC fc 

end for 

Allocation step: 

for all Xsi (t) with « = 1, . . . , n do 

find the cluster index k, for h* G [m k ; M k \: 

XsM -X C fe </ E,<^(7 s< fe W - 72C0) 2 P* < Efc< fc .(7?W - 7*< V VA; ^ fc' 
end for 



3 Dealing with simulated and real data 

The performance of the proposed clustering method was evaluated by studying the results obtained 
through the application on simulated and real datasets. 

3.1 Test on Simulated data 

First datasets are generated from a spatio functional random field with different spatial functional vari- 
ability structure. 

Specifically, given a sample of curves (x si (i), • • • , XsAt), ■ ■ ■ , Xs„(0) for t E T where s, is a generic 
data location in the rf-dimensional Euclidean space, and % s . (t) is generated by a spatio-functional Gaus- 
sian random field. 

The primary scope is to test the performances of the procedure in detecting spatio functional vari- 
ability structures. Thus, it is considered a situation largely used in geostatistics, where the covariance 
between Xsi(t), Xsj (t) is a stationary separable function of the form: 

Csep (M) = cov { Xsi (t),XsM} = °s (h) C T (u) (13) 

where C s (h) and Ct (u) are stationary, purely spatial and purely temporal covariance functions, respec- 
tively, defined on two generic locations Sj, Sj that are apart by h = Sj — Sj with a time span u — \U — t 3 -\. 

10 



The simulation schema proposed by (Sun, Genton (2011)) is considered as reference. In particular 
the spatial covariance function has the following form: 

C s {h) = {l-v) exp (-c \h\) + u5 h=0 (14) 

where c > controls the spatial correlation intensity, and v E (0, 1] is the nugget effect; the temporal 
covariance function is of the Cauchy type having the following form: 

C T {u) = {u + a^y 1 (15) 

where a E (0, 1] controls the strength of the temporal correlation and a > is the scale parameter in 
time. 

Six datasets made by n = 300 curves located on a regularly spaced grid have been generated. The 
following model is used: 

Xs {t) = n a (t) + e.(t) t E T (16) 

with mean /U s (t) = and e s (t) is a Gaussian random field with zero mean and covariance function as 
defined above. Each simulated dataset is made by curves belonging to three clusters C\, C 2 , C 3 . Each 
cluster includes 100 spatially adjacent curves generated according to the parameter sets in table[T| 

In each dataset and in each cluster there is no nugget effect (v = 0); moreover, the other parameters 
are set to a = 1 and a = 0, 1. 

There are two basic scenarios which are different in the values of standard deviation a used for 
generating the Gaussian random field of a cluster, so that the datasets 1, 2, 3 belong to the first scenario, 
while the datasets 4,5,6 belong to the second one. 

The datasets of both scenarios are designed to get three different levels of spatial correlation intensity 

c. 

In order to evaluate the capability of the proposed method to discover the spatial variability structures 
in the data and the curves which concur to form them, the well known Rand Index (Rand (1971)) is used. 
This index, whose value is in the range [0,1], allows the measurement of the degree of consensus between 
two partitions so that the value indicates that the two partitions do not agree on any pair of items while 
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Values of a 


Values of c 


Dataset Id 


Ci 


c 2 


c 3 


Ci 




c 3 


1 


5 


10 


15 


3 


7 


10 


2 


5 


10 


15 


5 


7 


9 


3 


5 


10 


15 


3 


9 


15 


4 


7 


10 


13 


3 


7 


10 


5 


7 


10 


13 


5 


7 


9 


6 


7 


10 


13 


3 


9 


15 



Table 1 : Parameters for simulated datasets 

1 means that the partitions are exactly the same. 

The test consists in computing the Rand Index between the true partition of data which emerges from 
the simulation schema and the partition given as output by the proposed clustering method. Since the 
latter depends on the initial random partitioning of data, the following table reports, for each dataset, the 
average Rand Index calculated on 100 repetitions of the algorithm. 



Dataset Id 


Average Rand Index 


1 


0.88 


2 


0.87 


3 


0.85 


4 


0.84 


5 


0.82 


6 


0.79 



Table 2: Rand Index value for each simulated dataset. 

The clustering results for the six datasets reflect the expectations based on the simulations. The RI 
appears to be high for all the simulated datasets, especially for the first dataset, where the value is 0.88. 
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Figure 1: Clustering results plotted on the spatial grid for the datasets 1,2,3. The color of the dots 
identifies the cluster membership. 

The results are very interesting, since the clustering structures in data are discovered. The good 
performance of the method is also highlighted by a graphic representation in Figure 1, 2, which plots 
the spatial locations of the three different clusters. Finally, Figure [3] highlights the different variability 
structures through clusters prototypes. 
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Figure 2: Clustering results plotted on the spatial grid for the datasets 4, 5, 6. The color of the dots 
identifies the cluster membership. 
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Figure 3: Theoretical variogram functions for the simulated datasets 

3.2 Test on real data 

In order to evaluate the performance of the proposed strategy on real data, a dataset was provided by the 
Institute for Mathematics Applied to Geo sciences^] The dataset reports the average monthly temperatures 
recorded by approximately 8000 stations located in the US, in the period 1895 to 1997. 

Tests used data from 1993 — 1997; thus for each station there is a time series made by a maximum 
of 60 observations. Since for several stations there are no data in the considered period, the dataset is 
composed of 4500 time series. 

2 http://www.image. ucar.edu/Data/US.monthly.met/ 
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The first step of the analysis is to construct the set of functions expanded in terms of 5-Spline Basis 
functions Q. An appropriate order of expansion Z is chosen, taking into account that a large Z causes 
overfitting and a too-small Z may cause important aspects of the function to be missing of the estimated 
function (Ramsay, Silverman (2005)). They consider a procedure based on a classical non-parametric 
cross-validation analysis. For each series, cubic splines are evaluated in order to produce a collection of 
smooth curves that is able to take into account the variability of the data. 

The very large extension of the spatial region involved in the monitoring activity makes it difficult to 
apply geostatistics methods based on the assumption of stationarity. Since stationarity and isotropy are 
assumed in the strategy the spatial trend is removed in a first step of the analysis by using a functional 
regression model with functional response (smoothed temperature curves) and two scalar covariates (lon- 
gitude and latitude coordinates in decimal degrees) (Giraldo et al. (2009)). 

On these spatially located curves, it is evaluated the capability of the proposed strategy in discovering 
different variability structures and their associated spatial regions. 

In order to run the clustering algorithm, the following input parameters have to be set: 

• the number of clusters K 

• the theoretical variogram model to fit the empirical one for each cluster 

Since there is not any information on the true number of spatial variability structures, the algorithm 
is applied for K = 2, . . . , 6 and then K is selected according to the maximum decreasing of the value of 
the optimized criterion A(P, L). For the tested dataset the best choice is K = 3. 

The theoretical variogram model is chosen evaluating several well known parametric models: Espo- 
nential, Spherical, Gaussian. The procedure is run for each model starting from the same initialization 
and then the fitting of each model to the data is evaluated, measuring the value of the criterion A(P, L). 
The results in Table [3] highlight that the best model is the exponential variogram thus, it is used on the 
tested dataset. 

Starting from the chosen input parameters, the algorithm run on the dataset, detects the spatial regions 
available in Fig. [4j The value of the optimized criterion is A(P, L) = 2.9e +4 ; the number of iterations 
until convergence is 9. 
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Trace-variogram 
model 
Exponential 
Gaussian 
Spherical 

Table 3: Criterion evaluation for several theoretical variogram models. 
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Figure 4: Clusters plotted on the geographical map. 

It is possible to note that the three discovered clusters split the studied area into three spatial regions, 
which include most of the east and west coasts, a northern area and a southern area. 

These spatial regions are characterized by three different spatial variability structures as shown in 
FigS 

1200 
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800 

400 
200 

°0 20 40 60 "0 20 40 60 "0 20 40 60 

Figure 5: Theoretical variogram models for the three discovered clusters. 

It is possible to note that the variogram corresponding to the third cluster shows the lowest level of 
variance (sill); the second cluster presents a variogram with highest sill level. The range of the variograms 
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A(P,L) 

2.9e+ 4 
3.5e+ 4 
3.6e+ 4 





is 29 for the first cluster, 25 for the second cluster and 16 for the third one. 

Looking at the plots, it is possible to note that the variability in the first and second clusters rises at a 
lower rate when it is compared to the third cluster. 

4 Summary and conclusions 

This paper has introduced an exploratory strategy for geostatistical functional data. 

It is a dynamic clustering method that partitions a set of geostatistical functional data into clusters 
that are homogeneuos in terms of spatial variability and that represents each cluster with a prototype 
variogram function. 

The approach is distinct from others since it discovers both the spatial partition of the data and the 
spatial variability structures representative of each cluster. The spatial information is incorporated into 
the clustering process by considering the variogram as a measure of spatial association, emphasizing the 
average spatial dependence among curves. 

This strategy can represent a very interesting methodological proposal for analyzing georeferenced 
curves in which spatial dependence plays an important role in exploring the similarity among curves. 
As in classical geostatistics data analysis, it assumes that the process generating data is stationary and 
isotropic. However, an alternative would be to consider an anisotropric process where the spatial depen- 
dence changes with the direction. In this case, it would be interesting to introduce a directional variogram 
model for functional data and demonstrate the main characteristics. 
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